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Abstract 

Liquid water, at ambient conditions, has short-range density correlations which are well known 
in literature. Surprisingly, large scale molecular dynamics simulations reveal an unusually long- 
distance correlation in 'longitudinal' part of dipole-dipole orientational correlations. It is non- 
vanishing even at 75 A and falls-off exponentially with a correlation length of about 24 A beyond 
solvation region. Numerical evidence suggests that the long range nature of dipole-dipole cor- 
relation is due to underlying fluctuating network of hydrogen-bonds in the liquid phase. This 
correlation is shown to give a shape dependant attraction between two hydrophobic surfaces at 
large distances of separation and the range of this attractive force is in agreement with experiments. 
In addition it is seen that quadrupolar fluctuations vanish within the first solvation peak (3 A) 

PACS numbers: 61.20.Ja, 61.25.Em, 82.30.Rs 



* Electronic address: 'jmpka nth@inisc.res.in| 
Electronic address: vani@inisc.res.inl 
•l-Electronic address: lramesha@imsc.res.inl 



1 



Water molecule with its hydrogens and lone pairs in tetrahedral arrangement makes 
hydrogen-bonds with its neighbouring molecules. In the liquid phase the hydrogen-bond 
pattern undergoes rapid fluctuations at pico-second time scales , U] , thus resulting in 

large orientational entropy. It is well known that this special property bestows liquid water 
with some unique properties, in particular the hydrophobic force of attraction between 
non-polar solutes. Clever experiments have been performed to measure quantitatively the 
distance properties of the hydrophobic force between mesoscopic surfaces js]. Understanding 
these distance properties is necessary initial step to develop a proper theory for bulk liquid 
water. Here we make a preliminary attempt towards the same using molecular dynamics 
(MD) simulations and general principles of statistical mechanics. 

A water molecule can be modelled as a set of five points corresponding to neutral oxygen 

0, two positively polarized hydrogens Hi, H2 and two negatively polarized lone-pair sites 
Li, L2 placed at tetrahedral angles about the oxygen atom. The angles between OH 12 and 
0Li^2 and the length of each of these vectors can fiuctuate. Such a molecule's orientations 
can be conveniently described with a choice of vectors defined as : 

^ OH1 + OH2 ,.0Li + 0L2 

''^'^^'^ = \0H, + 0H2\ - ^+^ |OL, + OL,| 

where 'r' is the position of oxygen atom in the bulk. The choice of ei(r) and e2(r) is 
such that they do not depend upon bond lengths of the molecule; they are symmetric with 
respect to hydrogens and lone-pairs of the molecule. ei(r), e2{r) and 63 = ei x 62 are the 
corresponding orthonormal vectors. Here ei(r) is dominantly the direction of dipole field 
and 62 (r) exists only if the water molecule differs from its mean (near-tetrahedral) geometry 

1. e. it is proportional to the quadrupole moment of the molecule. 

The e-vectors [Eq. ([1])] form a complete triad with which orientation of any vector {OH 
or OL ) can be specified. Consequently dynamics of water can be understood to be an 
interacting system of the e-vector fields. In particular MD simulation of water molecules 
implicitly gives us the dynamics of these fields. There upon various statistical correlations 
involving ei(r), e2{r) and p(r) = (ei(r))^ = (e2(r))^ in the liquid phase of water can be 
formulated as follows : 
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< p(ri)p(r2) > = g{r) (2a) 



r 



< p{ri)ea{r2) > = -da{r) (2b) 



r 

2^ r 



{6^^ - 3— )U(r) (2c) 



where r = (rl — r2), r = |f|; subscripts a,b = 1,2 (denote either of 61,62) and vector 
indices i,j = 1,2,3 (denote directions in three dimensional space). 

The translational and rotational symmetry of the system enables decomposing the ten- 
sorial properties of these correlations explicitly and thus analyze the data in terms of simple 
scalar functions like g{r), da{r), tab{r), lab{r). The function g{r) is the radial distribution 
function and it portrays distance-dependant density correlations only (here, of oxygen). The 
remaining functions capture the correlations among other degrees of freedom of the vector 
fields. 

TIP5P model possesses all orientational degrees of freedom of a water molecule and has 
improved accuracy in predicting the structural properties of water at ambient conditions. 
The simulations of TIP5P water system are performed with GROMACS (version 3.3.1) 
package 6| with an integration time step of 2 fs. The fast-moving bonds O — H are con- 
strained using LINGS algorithm. A large system consisting of 110592 molecules in a 150 
A box is equilibrated for 2 ns in constant pressure (isotropic and 1 atm) and temperature 
(300 K) NPT ensemble followed by a production run of 2 ns in a constant volume NVT 
ensemble. The configurations are saved every 100 ps for analysis. A cut-off distance of 
12 A and a pair-list distance of 15 A are used to compute all non-bonded interactions and 
periodic boundary conditions are imposed. Full electrostatic interactions are computed with 
Particle Mesh Ewald (PME) method with a tolerance of 10~^ and updated every two time 
steps 0, 8|. 

Density correlation g{r) of TIP5P displays all the well-known solvation peaks; in addition, 
due to large system size and hence 1 setter statistics, few more prominent peaks are observed 
at about r = 8.8 A and r = 10.8 A 9]. The orientations of (dipolar) field ei are analyzed by 
the correlations < ei*(0)ei-'(r) > [Eq. (j2])] where i,j refer to components of Ci vector. This 
is conveniently decomposed into two parts : transverse trace part tii(r) =< ei(0) ■ ei(r) > 
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which measures the dipoles' ahgnment with respect to each other and thus solely contributes 
to Kirkwood dielectric function fiol. [ill. 12|: and longitudinal traceless part /ii(r) =< ei(0) • 
f ei{r) ■ f > which is a measure of alignment of the vectors with respect to radial vector 
separating them. 

The transverse correlation function tii(r) shows oscillatory solvation structure, but van- 
ishes beyond 14 A [Fig. (Q] (in compliance with the rotational symmetry in the full system). 
The function In (r) is seen to be always positive and furthermore, in the 14 — 75 A regime 
it can be fitted to Ornstein-Zernike (OZ) form as 



g-r/5.2(l) e-^/24(i: 
/^^(r) = 0.39(2) + 0.027(1) - 



r 



r > UA (3) 

/ii(r) shows longest correlation length of 24 A. Furthermore it exhibits solvation peaks 
upto 14 A [Fig. (j2])]. In our simulation data upto 75 A the statistical sampling errors 
dramatically reduce as we go to large distances (as expected) j^. 

The dipole-oxygen correlation di{r) =< ei(0) ■ fp{r) > also exhibits solvation structure 
and vanishes beyond 14 A. It is also found that correlations involving 62, 63 all vanish upto 
statistical errors beyond the first solvation peak 9(|. Therefore 62, the quadrupole moment 
of water, fluctuates locally and randomly without any non-local correlations. 

TIP3P model [l^, by design, has ei degree of freedom only i.e. each water molecule's 
orientation can be completely described by ei field alone. The simulations on TIP3P wa- 
ter system are performed using NAMD (version 2.6) [14]. Here, 33105 water molecules are 
simulated in a cubical box of size 100 A and the procedures employed for collecting equilib- 
riated configurations are same as those described in case of TIP5P. The constrained model is 
implemented using SETTLE algorithm. Analysis in this case too shows that tii(r) vanishes 
beyond solvation region, whereas lu (r) follows the same asymptotic behaviour as described 
by Eq. m- 

A water molecule in liquid phase is predominantly influenced by hydrogen-bonding (short- 
range interaction) and further, it has a net dipole moment, which interacts through long- 
range Coloumbic forces. We would like to ascertain if the long distance behaviour of lii{r) is 
due to the short-range hydrogen-bond interactions or the long-range Coloumbic interactions 



15|, 



16]. To test this possibility, the Coulombic interactions are smoothly truncated at 12 A 



in TIP3P model, thus retaining an effective short-range interaction alone which imitates the 
effect of hydrogen-bonding. We find that lu (r) remains essentially unchanged in the regions 
of first few solvation shells and r > 30 A. The intermediate region exhibits over-structuring 
effects upto 30 A 0, 3 • 

The above three cases are in agreement with Eq. asymptotically. These observations 
suggest that (i) water in liquid phase has fluctuations only in dipole degree of freedom; in 
contrast the quadrupole has no effect beyond the first solvation peak, (ii) these dipole fluc- 
tuations in liquid water are influenced by local environment of respective molecule, through 
hydrogen-bonding, significantly more compared to long-range electrostatic interactions, (iii) 
furthermore, the dipole fluctuations exhibit long distance correlations. Simulations show 
that OZ like behaviour is exhibited only by longitudinal dipole-dipole correlation and not 
transverse component or oxygen-oxygen density correlation. Recently there has been an 
experimental observation which reported correlation length of the order of a nanometer in 
liquid water 191]; in line with the shorter correlation length in /ii(r). Below we suggest 
that the longer correlation length of 24 A is indirectly observed in SFA experiments as force 
between hydrophobic plates. 

Hydrophobic effect : The first notable mechanism postulated to describe the origin of 
hydrophobic force came from solvation studies of Frank and Evans in the name of "iceberg" 
model {20! and later, the same effect has been elucidated by Kauzmann on its possible biolog- 
ical implications [2l| . This phenomenon, in its various manifestations, has been extensively 



discussed in recent literature 
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251]. Experimentally the force of attraction between 



two nominally hydrophobic macroscopic surfaces has been measured and it was found that 
in the range 10 A to 100 A the force falls-off exponentially with a correlation length of 12 
A 26]. Later there have been several such studies using Surface Force Apparatus (SFA) 



with surfaces prepared and characterized using wide range of techniques [5 



were very few theoretical developments (Lum et al 



and ref 's therein, 129 



. Yet there 
30|) to explain 



qualitatively different force profiles observed in experiments. We address here monotonic 
nature of the force as observed in several experiments [s], Q, 27]. 

Hydrophobic surfaces cannot form hydrogen-bonds with the surrounding water, conse- 
quently water molecules rearrange themselves such that they form a sheet of hydrogen-bond 
network on the surface. Their interactions are such that the directions of lone pairs and 
hydrogen atoms are perpendicular to the surface normal of the hydrophobe. Owing to the 
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approximate tetrahedral conformation, water molecules cannot have a unique configuration 
satisfying the above criterion IsiJ. Consequently the y ex plore other possible orientations 
as well by fluctuating at the pico-second time scales 32, l33|]. These network fluctuations 
contribute significantly to the free energy of solvation of the hydrophobe. In presence of two 
such hydrophobes, as noted earlier, the range of the force acting between them is large and 
due to limited computational resources it is not possible to directly simulate and observe 
this effect numerically. Alternatively a quantitative theoretical estimate is being considered 
below. 

Interaction between hydrophobic surface and solvent water can be written in terms of 
n(r), the local unit normal vector to the hardcore van der Waals surface of the hydrophobe 
and ei(r'), the dipole of water molecule near the surface, where r' = r + 5r; 5r is typical 
length of hydrogen arm of water molecule (about 1 A). A simple local interaction term can 
be taken as (n(r) ■ ei(r'))^ implyin g th at the water dipoles orient orthogonal to the surface 



normal as seen in simulations |34l . l35l . l36| (importantly, no linear term in n ■ ei, for that 
means a preferential orientation of the water dipole inward/outward to the surface). 

The change in free energy due to purely hydrophobic interaction between two small 
surfaces 5*1 and 5*2 [Fig. ([3])] in water can be estimated by. 



AH = ^ dMMri) ■ ei{r[))^ 
^ Jsi 

+ ? / dn2{h2{r2) ■ Ciir'^))^ 

^ JS2 



e 



-AG/feT ^ ^ ^-AH/kT ^ 



where 7 is a measure of strength of interaction between hydrophobic solute and water 
which can depend upon temperature, density and other parameters defining the thermody- 
namic system. The brackets < ... > refer to statistical averaging with respect to pure water 
system and integration is over area of each surface. As illustrated in Fig. Q, Si, S2 refer to 
two arbitrary hydrophobic surfaces and i? is a vector along minimum distance of separation 
between them. 

When the distance R (= \R\) is large compared to radius of curvature of each surface 
and the surface areas sufficiently small, the statistical averaging can be done by cumulant 
expansion. The leading term of the force F{R) = —dAG/dR is given by the following 
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equation [9| 



< [ni(ri) ■ ei(r'i) n2{r2) ■ ei^ri^)]'^ >] 
= Wt^'"^' ^Tr[Ss,E(i?)Ss,E(i?)] (5) 
where Ai, A2 are areas of the surfaces and the matrices E(i?), Ss are given by, 



- -^(^'' - 3^)/n(i?) for large R 

= ^ [ dn n'n^ 
A Js 

For a segment of spherical surface (such as a hardcore van der Waals surface) subtending 
a cone angle 6 at its centre and A^* = M'/\M\ where = \ f^dh n\ t!^ = l6'^ - 1(6'^ - 
3N'N^)cos9{l + cos9). For a hemi-sphere, 9 = ^, T.'^ = ^6'^ . For a plane, ^ = 0, S'^ = N'N^ 

The above result on hydrophobic force is true very generally. As discussed in earlier 
paragraphs, the leading order (n-ei)^ is taken to be the interaction energy term for simplicity. 
By including the non-leading terms in the interaction energy function [Eq. (jl])] and doing 
the cumulant expansion, it can be shown that the force term [Eq. (^] for large R remains 
unchanged, thus establishing the generality of the result. 

These considerations are valid for distances beyond the solvation region of a typical water 
molecule. The cumulant expansion allowed decomposing the force equation as a simple 
convolution of surface-dependant part and water-dependant part. Eq. enables us to 
conclude that range of the force between hydrophobic surfaces at large distances is always 
attractive governed by l1i{R) — e'^^^"^ for large R. Therefore the hydrophobic force falls off 
exponentially with a largest correlation length of about 12 A, in addition to several other 
shorter range exponents as well. 

F{R) ~ -e-^/^'^ for large R (6) 

The strength of attraction is proportional to area A and shape of each surface given by 
the tensor S, the second moment of surface normal. The final trace operation over the 
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matrices E(R) and Ss [as in Eq. ([5])] implies that the hydrophobic attraction is not just 
a purely distance dependant interaction like van der Waals'. Indeed the orientation of the 
surface shapes relative to each other can modify this force significantly. As an example if 
two small planar hydrophobic surfaces are mutually perpendicular and are sufficiently far 
apart, there should be no force between them as opposed to when they face each other. 



[1] F. H. Stillinger, Science 209, 451 (1980). 

[2] F. W. Starr, J. K. Nielsen, and H. E. Stanley, Phys. Rev. Lett. 82, 2294 (1999). 
[3] F. N. Keutsch and R. J. Saykally, Proc. Natl. Acad. Sci. U. S. A. 98, 10533 (2001). 
[4] C. J. Fecko et al, Science 301, 1698 (2003). 

[5] E. E. Meyer, K. J. Rosenberg, and J. N. Israelachvili, Proc. Natl. Acad. Sci. U. S. A. 103, 
15739 (2006). 

[6] E. Lindahl, B. Hess, and D. van der Spoel, J. Mol. Mod. 7, 306 (2001). 

[7] D. Frenkel and B. Smit, Understanding Molecular Simulation ( Computational Science Series, 

Vol 1) (Academic Press, U.S.A., 2001). 
[8] A. Leach, Molecular Modelling: Principles and Applications (2nd Edition) (Prentice Hall, 

U.S.A., 2001). 

[9] See supporting information at www. imsc. res. in/'^maruthi/research. html 
[10] J. G. Kirkwood, J. Chem. Phys. 7, 911 (1939). 

[11] P. L. Silvestrelh and M. Parrinello, Phys. Rev. Lett. 82, 3308 (1999). 
[12] M. Sharma, R. Resta, and R. Car, Phys. Rev. Lett. 98, 247401 (2007). 
[13] W. L. Jorgensen et a/., J. Chem. Phys. 79, 926 (1983). 
[14] L. Kale et al, J. Comput. Phys. 151, 283 (1999). 

[15] J. -P. Hansen and I. R. McDonald, Theory of Simple Liquids, Third Edition (Academic Press, 
U.K., 2006). 

[16] G. Mathias and P. Tavan, J. Chem. Phys. 120, 4393 (2004). 

[17] M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids, Oxford Science Publications, 

1 ed. (Clarendon Press, Oxford, 1987). 
[18] J. Kolafa and I. Nezbeda, Mol Phys 98, 1505 (2000). 

[19] C. Huang et al, Proc. Natl. Acad. Sci. U. S. A. Epub before print doi: 



8 



10.1073/pnas.0904743106 (2009), (The authors attribute the correlation length to density 

fluctuations. This aspect needs to be investigated in more detail.). 
[20] H. S. Prank and M. W. Evans, J. Chem. Phys. 13, 507 (1945). 
[21] W. Kauzmann, Adv. Protein Chem. 14, 1 (1959). 
[22] P. Ball, Chem. Rev. 108, 74 (2008). 
[23] M. Chaplin, Nat. Rev. Mol. Cell Biol. 7, 861 (2006). 

[24] K. A. Dill, T. M. Truskett, V. Vlachy, and B. Hribar-Lee, Annu. Rev. Biophys. Biomol. Struct. 

34, 173 (2005). 
[25] D. Chandler, Nature 437, 640 (2005). 
[26] J. N. IsraelachviH and R. Pashley, Nature 300, 341 (1982). 

[27] H. K. Christenson and P. M. Claesson, Adv. Colloid Interface Sci. 91, 391 (2001). 

[28] K. Lum, D. Chandler, and J. D. Weeks, J. Phys. Chem. B 103, 4570 (1999). 

[29] F. Despa, A. Fernandez, and R. S. Berry, Phys. Rev. Lett. 93, 228104 (2004). 

[30] N. A. M. Besseling, Langmuir 13, 2113 (1997). 

[31] Y. K. Cheng and P. J. Rossky, Nature 392, 696 (1998). 

[32] Y. L. A. Rezus and H. J. Bakker, Phys. Rev. Lett. 99, 148301 (2007). 

[33] A. Poynor et al, Phys. Rev. Lett. 97, 266101 (2006). 

[34] P. L. Chau, Mol. Phys. 99, 1289 (2001). 

[35] P. Jedlovszky, J. Phys.: Condens. Matter 16, 5389 (2004). 

[36] T. M. Raschke and M. Levitt, Proc. Natl. Acad. Sci. U. S. A. 102, 6777 (2005). 



9 




10 



0.0075 



0.005 



0.45 
0.3 
0.15 






10 



15 



60 



FIG. 2: Exponential decay in longitudinal dipole-dipole correlation of liquid water outside 

the solvation region, {lower) TIP5P data and fit function, given by Eq. ([3|), right on top of each 
other, (middle) TIP3P data, (upper) TIP3P with truncated Coulombic interactions. For clarity, 
the middle and upper plots are shifted up by 0.001 and 0.002 units respectively, (inset) Zii(r) 
inside the solvation region. 




FIG. 3: 5i, 5*2 are hydrophobic surfaces with their local normal vectors ni, h2- '-R' is the minimum 
distance between the two surfaces 
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